Testing LIGER

# library(devtools)
# install_github('MacoskoLab/liger')
library(tidyverse)

── Attaching packages ───────────────────────────────────────────────────────────────────────────────────── tidyverse 1.2.1 ──
✔ ggplot2 3.2.1     ✔ purrr   0.3.2
✔ tibble  2.1.3     ✔ dplyr   0.8.3
✔ tidyr   1.0.0     ✔ stringr 1.4.0
✔ readr   1.3.1     ✔ forcats 0.4.0
── Conflicts ──────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ tidyr::expand() masks Matrix::expand()
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
✖ tidyr::pack()   masks Matrix::pack()
✖ tidyr::unpack() masks Matrix::unpack()

library(liger)
library(Rtsne)
library(reshape2)

Attaching package: ‘reshape2’

The following object is masked from ‘package:tidyr’:

    smiths

library(ggrepel)

Testing LIGER for integration of scATAC and scRNA data. Following vignette from (gitHub)[https://macoskolab.github.io/liger/walkthrough_rna_atac.html].

1. Download and load datasets:

ATAC data is preprocessed to gene-level features (no. of reads that overlap gene body and promoter for each gene). Data is downloaded from (here)[https://umich.app.box.com/s/5hkhpou4inrulo570yhc38ay8g3ont5b]

rna_clusts = readRDS("~/10X_data/rna_cluster_assignments.RDS")
atac_clusts = readRDS("~/10X_data/atac_cluster_assignments.RDS")
pbmc.atac <- readRDS('~/10X_data/pbmc.atac.expression.mat.RDS')
pbmc.rna <- readRDS('~/10X_data/pbmc.rna.expression.mat.RDS')

2. Create LIGER object

List of count matrices as input. createLiger converts to sparseMatrices.

pbmc.data = list(atac=pbmc.atac[,names(atac_clusts)], rna=pbmc.rna[,names(rna_clusts)])
int.pbmc <- createLiger(pbmc.data)

[1] "Removing 97 genes not expressing in atac."
[1] "Removing 11048 genes not expressing in rna."

3. Data preprocessing

## Normalize to column total (adj. for sequencing depth)
int.pbmc <- normalize(int.pbmc)
## Select highly variable genes ONLY IN THE RNA DATASET
int.pbmc <- selectGenes(int.pbmc, datasets.use = 2)
## Scale to the SD
int.pbmc <- scaleNotCenter(int.pbmc)

4. Run NMF

int.pbmc <- optimizeALS(int.pbmc, k=20)

  |                                                                                                                           
  |                                                                                                                     |   0%
  |                                                                                                                           
  |=                                                                                                                    |   1%
  |                                                                                                                           
  |==                                                                                                                   |   2%
  |                                                                                                                           
  |====                                                                                                                 |   3%
  |                                                                                                                           
  |=====                                                                                                                |   4%
  |                                                                                                                           
  |======                                                                                                               |   5%
  |                                                                                                                           
  |=======                                                                                                              |   6%
  |                                                                                                                           
  |========                                                                                                             |   7%
  |                                                                                                                           
  |=========                                                                                                            |   8%
  |                                                                                                                           
  |===========                                                                                                          |   9%
  |                                                                                                                           
  |============                                                                                                         |  10%
  |                                                                                                                           
  |=============                                                                                                        |  11%
  |                                                                                                                           
  |==============                                                                                                       |  12%
  |                                                                                                                           
  |===============                                                                                                      |  13%
  |                                                                                                                           
  |================                                                                                                     |  14%
  |                                                                                                                           
  |==================                                                                                                   |  15%
  |                                                                                                                           
  |===================                                                                                                  |  16%
  |                                                                                                                           
  |=====================================================================================================================| 100%
Converged in 3.941199 mins, 16 iterations.
Best results with seed 1.

smp <- sample(1:nrow(int.pbmc@H$atac), size = 1000)
smp_genes <- sample(1:length(int.pbmc@var.genes), size = 100)
int.pbmc@W[,smp_genes] %>% heatmap(Rowv = NA)

Error in int.pbmc@W[, smp_genes] %>% heatmap(Rowv = NA) : 
  could not find function "%>%"

Align factor

Dodgy quantile normalization

int.pbmc <- readRDS(file = "~/10X_data/pbmc.trainedLIGER.RDS")

cannot open compressed file '/home/jovyan/10X_data/pbmc.trainedLIGER.RDS', probable reason 'No such file or directory'Error in gzfile(file, "rb") : cannot open the connection

int.pbmc@scale.data$atac[,
                         smp_genes] %>%
  melt(varnames=c("cell", "gene")) %>%
  ggplot(aes(value, color=gene)) + geom_histogram() +
  facet_wrap(gene~.)
int.pbmc@scale.data$rna[,
                         smp_genes] %>%
  melt(varnames=c("cell", "gene")) %>%
  ggplot(aes(value, color=gene)) + geom_histogram() +
  facet_wrap(gene~.)

  

Visualize dataset-specific VS common latent factors

Even in dataset specific component there is separation of cell types (Meaningful information that is ATAC specific?)

atac.comp <- int.pbmc@H$atac[smp,] %*% int.pbmc@V$atac

Error: object 'smp' not found

rna.comp <- int.pbmc@H$rna[smp,] %*% int.pbmc@V$rna
rna.common.comp <- int.pbmc@H$rna[smp,] %*% int.pbmc@W
tsne.rna.comp <- Rtsne(rna.comp, pca=F, check_duplicates = F, theta=0.5, perplexity=30)
rownames(tsne.rna.comp$Y) <- rownames(rna.comp)
tsne.rna.common.comp <- Rtsne(rna.common.comp, pca=F, check_duplicates = F, theta=0.5, perplexity=30)
rownames(tsne.rna.common.comp$Y) <- rownames(rna.common.comp)
tsne.rna.comp$Y %>%
  as.tibble(rownames="cell") %>%
  mutate(clust=rna_clusts[cell]) %>%
  ggplot(aes(V1, V2, color=clust)) + 
  geom_point()

tsne.rna.common.comp$Y %>%
  as.tibble(rownames="cell") %>%
  mutate(clust=rna_clusts[cell]) %>%
  ggplot(aes(V1, V2, color=clust)) + 
  geom_point()

Compare weights of dataset specific and common factors

Are HVGs in RNA low coverage in ATAC?

Evaluate feature selection for ATAC-seq data